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We consider a rotating Bose-Einstein condensate in an optical lattice in the regime in which the 
system Hamiltonian can be mapped onto a Josephson junction array. In an approximate scheme 
where the couplings are assumed uniform, the ground state energy is formulated in terms of vortex 
configuration. Application of method for ladder case presented and the results are compared with 
Monte-Carlo method. 
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I. INTRODUCTION 



II. ROTATING BEC IN AN OPTICAL LATTICE 



After the first experimental realization of the Bose- 
Einstein Condensates (BECs) 1 -, this field and its related 
topics attracted more attentions^. BEC in an opti- 
cal lattice 4 ^ has relation with many problems in the 
condensed matter physics, e.g. Bloch Oscillations^— , 
Wannier-Stark Ladders 1 ^, Josephson junction arrays 1 - 1 - - — 
and superfluid to the Mott-insulator transition^ 8 . - — . 

Behavior of a rotating BEC in an optical lattice, 
is similar to that of a superconductor in the mag- 
netic fielder—. Here with a trap potential without the 
lattice structure, the Abrikosov vortex lattice can be 
observed 1 ^ 1 ^ ' 14 ' 15 £L2L , When the lattice structure adds 
to the trap potential, the vortex lattice changes 3 ^ and 
with some criterions, this study can be mapped on the 
problem of the Josephson junction arrays (JJAs); there 
will be the same vortex structure for both system o 12 ' 14 . 
Most of the studies in this field has been done on the 
lattice potentials with the square symmetries and few 
works on the BEC with other symmetric a 33 ' 34 . Yet due 
to the possibility of the realization of the optical lattices 
with different symmetries^and with more than one spa- 
tial frequency 3 ^, the study of the BEC in the quasiperi- 
odic potentials can be usefu l 33 ' 34 . 



In this section we show how in some approximation, 
rotating BEC in an optical lattice can be modeled as 
Josephson junction array a 11 ' 12 ' 14 . Hamiltonian for a BEC 
in an optical lattice with a rotation frequency flz is, 



H = / dr¥[ 



(—iKV — mQz x r) 
2m 



V R , 



(2.1) 

where m is the atomic mass and g = 4irH 2 a/m the 
coupling constant with the s-wave scattering length a. 
Conservation of the total number of particles is en- 
sured by chemical potential ji. The external potential 
Vext — Vh + Vo consists of two parts: modified harmonic 
potential Vh = m(w 2 - £l 2 )r 2 /2 + mcj 2 z 2 /2 and the po- 
tential of the optical lattice Vo , which may be chosen as 
periodic* 3 ^, quasiperiodio 33 ' 34 , or random 4 ^. For exam- 
ple, the potential with square symmetry can be written 
as Vo = Vo[sin 2 (kx) + sin 2 (ky)], with the periodicity 

7T /k. 

For a large lo z , we can suppose that system is frozen in 
axial direction. If energy due to interaction and rotation 
is small compared to the energy separation between the 
lowest and first excited band, the particles are confined 
to the lowest Wannier orbitals 1 ^. Therefore We can write 
ip in the Wannier basis as, 



Here we focus on the study of the problem in the JJA 
regime. We formulate the problem of the ground state for 
the JJA Hamiltonian in terms of the vortex configuration 
and then we discuss about the exact numerical result for 
the few number of the lattice points, and the Monte- 
Carlo method for a large lattice. We exploit harmonic 
approximation for the cosine Hamiltonian of the system 
and find energy for any given vortex configuration. Also 
we study the problem for the ladder case; this case can 
be used as the reduced two dimensional case 3 - 7 - and can 
inspire the general behavior of the ground state structure 
in the 2D lattice case. This method is faster than the 
Monte-Carlo method for the original Hamiltonian and 
directly results the vortex lattice, while in other methods 
instead of the vortex, the circulation has been used for 
the determination of the vortex position 38 . 
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A(r').dr',) (2.2) 



where A = flz x r is the analog of the magnetic vec- 
tor potential, Wj(r) = w(r — r^) the normalized Wannier 

wave-function localized in the i-th well and hi — a\a,i the 
number operator. Substituting 12.21 in Hamiltonian 12.11 
leads to the Bose-Hubbard model in the rotating frame 



ff = -i V (a\aje~ Al3 + h.c.) 



<»j> 
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^(e* - n)hi + — ^hi{hi - 1), (2.3) 
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where means i and j are nearest neighbors, t ss 

— J drw*(— h 2 V 2 /2m + Vo)wj is the hopping matrix el- 
ement, tj = J drw*VhWi an energy offset of each lattice 
site and U — g J dr|wj| 4 the on-site energy. The rota- 
tion effect is described by Ay proportional to the line 
integral of A between the i-th and j-th sites: Ay = 
(m/H) Q A[r').dr'. 

If the number of atoms per sites is large (n, 3> 1), 
the operator can be expressed in terms of its amplitude 
and phase, and the amplitude by the c-number as di ~ 
yfnle %6i . Using hi = rij — id/dQi and 9i = 9 i: Eq. 12.31 
reduces to the quantum phase model 

h = - E JiicoWi - h + A ») - \ E w> {2A) 

<i,j> i 1 

where Jy = 2t^/rTih~j. Also, the atoms are assumed 
to be distributed such as they satisfy the condition 
£j + Urii — fi. The magnitude of Jy decreases from the 
central sites outward because n.; has the profile of an in- 
verted parabola as we can see from solution of ID case 
in Fig. [U 
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FIG. 1: Solution of GP equation for a periodic one dimen- 
sional potential. It can be seen that |^| 2 is enveloped by an 
inverted parabola shape. 

Eq. 2.4 is just the Hamiltonian of a lattice of small 
Josephson junctions with inhomogeneous coupling con- 
stants Jy. With Jy ^> U, we can neglect the kinetic 
energy term; Moreover, we suppose that coupling con- 
stants are equal Jy = J. This regime can be achieved 
when the depth of optical lattice Vb is between 18 and 
25.Br (Eh = ti 2 k 2 /2m is recoil energy) and average par- 
ticle number is h « YIQ^L. Then we have: 

H = - E cos ^l (2-5) 

<i,j> 

with 7y = 9i — 9j — Ay ■ Our purpose is to find minimum 
of this Hamiltonian as a function of condensate phases 
9i with two constraints as follows. First constraint origi- 
nates in the single valuedness of the wave function which 
leads to the quantization of vorticity: 

E 7« = 27r(n fc - a k f) (2.6) 

where sum is over fcth plaquette, a k is area of fcth plaque- 
tte, and n& is integer. In the case of a condensate with the 



rigid-body rotation, the frustration parameter is defined 
as / = 2f2 / k with the quantum circulation k — h/m. In 
two dimensional arrays, frustration represents density of 
vortices per number of lattice plaquettes 42 . 

The second constraint comes from conservation of the 
particles numbers. The particle current from the nodes 
can be assumed as the derivation of the energy function 
with respect to the corresponding phased. This leads to 

E»»(7«)=0, (2.7) 

where sum is over all j which is connected to ith node. 

Since Hamiltonian is periodic, can be [ a fe/J or 
[a/c/J + 1, and we can replace them by which is zero 
or one, and replace a k f with akf — L°fc/J • These two set 
of equations can completely determine all 7y for a given 
set of rc.fc. Also, for a lattice, number of plaquettes plus 
number of nodes is equal to number of connections minus 
one£% therefore one of the equations (e.g. one of particle 
conservation equations) can be neglected. Assuming the 
particle flow from the nodes are small, we can linearize 
the Eq |2.7l i.e. siwyij ~ 7y. Now, we have a linear set 
of equations which must be solved for a given set of in- 
tegers rifc, i.e. for a vortex configuration. Energy of this 
vortex configuration is denoted by E{ n y. Then instead of 
original minimization scheme, we can use the set of these 
energies for minimization of Hamiltonian. It means that 
minimum of Hamiltonian correspond to a vortex configu- 
ration which minimizes these set of energies, with number 
of vortex configuration equal to 2 N where N is number 
of plaquettes in the lattice. 



III. LADDER CASE 

In this section, we apply the formulation of previous 
section on a ladder geometry and find an exact analytical 
formula for energy of vortex configuration. It is known 
that a ladder geometry can give an insight about ground 
state vortex structure in 2D lattices^. As we will see 
the symmetry of the problem in this case, reduces the 
number of equations such that there will be N linear al- 
gebraic equations to solve, where N is number of plaque- 
ttes. One of difference between problems in Josephson 
junction arrays and BEC lattices is boundary condition 
which is imposed on lattice. In the Josephson junction 
arrays, the behavior of the arrays with large number of 
plaquettes is interesting 4 ^, therefore the boundary condi- 
tion of lattice is usually assumed as periodic, i.e. the first 
plaquette is connected to the last plaquette. For BECs, 
the number of lattice sites is usually finite and free or 
fixed boundary condition are more suitable. 

Conservation of current 12.71 imposes that current in 
upper and lower junctions of a plaquette and therefore 
their corresponding 7y will be equal. We denote gauge 
invariant phase difference for upper and lower junctions 
of fcth plaquette by 7^ and those of the vertical junctions 
between k and fc + 1th plaquette by jk.k+i- Then the 
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node current equations read 

sm(7fe,fe+i) = sinfrk) - sin(jk+i), (3-1) 
and the linearized approximation gives 



7fe,fe+l = Ik - lk+l- 



(3.2) 



Imposing this equation on the flux quantization 2.6, 
we have 



47fe - 7fc+i - 7fc-i = 2n(n k + [a k f\ - a k f), (3.3) 

for the fcth plaquette with fixed boundary condition: 
7jv+i = 7o = 0, where N is number of plaquettes in 
ladder. We have left now with TV equations and 2N vari- 
ables. Denoting coefficient matrix by C, we can write 



7fe 



N 



k=l 



2n(N k - A k f), 



(3.4) 



where bi = 2it(ni — aif), and C 1 is inverse of coefficient 
matrix and its elements are, 



^ki 



iV-|fc-I|+l 



■i 



\k-X\-N-X _ £N+l-k-l 



k+l-N-l 



_ £—l}(£N+l _ £-iV— 1 



(3.5) 

where £ = 2 + and A k and N k are defined with 
appropriate summation. 

By using Eq. 13.41 we can find the energy (normalized 
by the number of junctions) for a vortex configuration 
{n} as: 

I N 

E {n} = 3jv T l \52{ 2cos (Hk) + cos(j k - 7fc+l)} 



fc=l 



cos(7at)}. 



(3.6) 



Related to the problem of minimizing the Hamilto- 
nian, each vortex configuration can be vortex configu- 
ration with minimum energy for a definite value, or an 
interval of /. Then for each value of frustration we need 
calculate the energy for all the vortex configurations. As 
an example, we begin with a square ladder (aj = 1) with 
8 plaquette in which we just need to study < / < 1. 
In Fig. [2] we have plotted ground state energy by the 
above procedure: / is incremented from to 1 by 1/128 
and for each value of /, energy of all the configurations is 
calculated and the least value has been chosen. We have 
specified each vortex configuration by a number m which 
in the binary representation with N digit (N is number 
of plaquettes) gives the configuration of the vortices with 
every one, meaning there is vortex on that plaquette. For 
example m = 3 means a vortex configuration 00000011, 
i.e. there are two vortices on the 7th and 8th plaquettes. 

Fixing / we can demonstrate the dependence of the 
energy to the configuration of the vortices. In Fig. [3] we 



have plot energy vs. vortex configuration for / = 1/2 
and 2/5 where the vortex configuration is labeled as we 
noted above. For / = 1/2 we see a mirror symmetry, 



it means that vortex configurations m and 2 



have 



the same energy. More general rule for square ladder with 
arbitrary / can be deduced easily from above equations: 
vortex configuration m for / and 2 N — m for 1 — / has 
same energy. Therefore, if vortex configuration m has the 
minimum energy for /, then vortex configuration 2 N — m 
has the minimum energy for 1 — /. This means that for 
/ = 1/2, the minimum energy is twofold degenerate. For 
example for the ladder considered above, both 01010101 
and 10101010 configurations give the minimum energy. 
In Fig. 2]we have shown position of vortices for different 
/. This figure shows that behavior of vortex is similar to 
the behavior of a charge q on a lattice with opposite sign 
onsite charge o 42 ' 44 . When we have one vortex on lattice, 
the vortex prefers to sit on middle of ladder, where in 
this case it can be two plaquettes 4 or 5. When we have 
two vortices on ladder, they prefer to divide ladder into 
three equal parts, which means that they must be on 
plaquettes 3 and 6. This description can be applied to 
vortex configurations with more vortices. 




FIG. 2: Ground state energy of a square lattice with 8 plaque- 
ttes. For each value of /, energy of all vortex configurations 
is calculated and the smallest value has been chosen. 
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FIG. 3: Energy of vortex configurations for square ladder with 
8 plaquettes. Binary representation is used to demonstrate 
vortex configurations as is explained in the text. / = 1/2 in 
(a) and / = 2/5 in (b). 

For more accurate results, we can relax linear approxi- 
mation for the current equation and use Newton- Raphson 
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FIG. 4: Position of vortex for a square ladder with 8 pla- 
quettes. Vertical axis shows the position of plaquette in the 
array. 



method, with the initial estimation from linear approx- 
imation. For above case, we compare linear approxima- 
tion and 10 iterates of Newton-Raphson method in Fig. 
[5l It worth to mention while energy differences are small 
but, in some small intervals of /, different vortex config- 
urations give the minimum energy. 



To check the validity of our method, instead, we apply 
the method on a quasiperiodic lattice e.g. a Fibonacci 
ladder. The Fibonacci ladders contains two types of 
plaquettes with lengths 1 and r = (1 + -\/(5))/2. If 
we denote plaquette with length 1 and cr, by S and 
L respectively, then these two types of plaquettes ar- 
ranged in bases of the following rule: for (n + l)th 
step of construction, we add (n — l)th step to end of 
nth step. With first step S and second step L, we 
have: LS, LSL, LSLLS, LSLLSLSL, .... The results arc 
shown in Fig. [5] along with the result of using the Monte- 
Carlo on the original Hamiltonian, with a good agree- 
ment (we note that for a aperiodic lattice it is not suffi- 
cient to study < / < 1). For a square lattice, increasing 
/ from to 1/2, both the number of the vortices and the 
magnitude of vorticity grow- 4 ^. and the linear approxima- 
tion may seem unreasonable. Fig. [6] shows that despite 
this fact, the approximate results are in a good agreement 
with the results of the direct method. 



IV. MC FORMULATION FOR VORTEX 
CONFIGURATION 

When N number of plaquettes grows, number of pos- 
sible vortex configurations grows as 2 , and expect for 
small number of plaquettes we can not determine the ex- 
act minimum vortex configuration which needs checking 
all of vortex configuration. Therefore we use a Monte- 
Carlo method: we start from a high temperature T and 
a random vortex lattice n and calculate its energy E n . 
Then we decrease temperature gradually and for each 
temperature we do following process a few times: we 
change the configuration randomly to n' with energy E n i , 
we accept the new configuration if exp(~AE/T) < r 
where r is a random number between and 1, and 
AE = E n — E n >. When temperature receives 0.0001, 
we take the final configuration as ground state. 




FIG. 5: Comparison between linear approximation (blue) and 
the Newton-Raphson method (green). 

For a periodic lattice finding the ground state from 
the MC method is troublesome for nonzero frustration 45 . 




FIG. 6: Comparison between energy obtained by the direct 
Monte-Carlo method (MC) and by the Monte-Carlo method 
with vortex configuration formulation (MCVOR). 



V. CONCLUSION 

The configuration of the vortex lattice of rotating 
BECs in an Optical lattice is investigated. For deep Op- 
tical Lattices when the number of particles in each site is 
large, the problem of rotating BEC in an optical lattice, 
can be mapped onto the model of arrays of Josephson 
junctions in presence of an external magnetic field. In 
this approximation we have formulated a vortex config- 
uration for the ground state. Our result for ladder case 
is presented and is in good agreement with Monte-Carlo 
result with original JJA Hamiltonian. This method can 
be extended to 2D case, and also cases with non- uniform 
coupling which looks more relevant for study of rotat- 
ing BECs. In these cases, with linear approximation, we 
deal with a set of linear algebraic equations for each vor- 
tex configuration which can be solved easily by use of 
proposed method. 

Effect of non-uniform coupling can affect the results 
especially in the situation of high density of vortices as 
in the case / = 1/2 for square lattice. Also as a better 
approximation we can consider finite U , then problem is 
similar to the arrays of small Josephson junctions 4 ^— . 
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We can also suppose that coupling is time-dependent 
which can occur in case of a vibrating optical lattice^. 



Then problem can be treated as a set of coupled pendu- 
lum equations with time-dependent lengths^. 
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